home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Turnbull China Bikeride
/
Turnbull China Bikeride - Disc 1.iso
/
ARGONET
/
PD
/
MATHS
/
RLAB
/
RLAB125.ZIP
/
!RLaB
/
toolbox
/
ode78
< prev
next >
Wrap
Text File
|
1994-07-18
|
4KB
|
160 lines
// ODE78 Integrates a system of ordinary differential equations using
// 7th order formulas.
//
// out = ode78(F, t0, tfinal, y0, Uoutput, tol, trace)
//
// INPUT:
// Ufun - variable containing user-supplied problem description.
// Call: yprime = fun(t,y)
// t - Time (scalar).
// y - Solution row-vector.
// yprime - Returned derivative vector; yprime(i) = dy(i)/dt.
// t0 - Initial value of t.
// tfinal- Final value of t.
// y0 - Initial value column-vector.
// Uoutput User output function
// tol - The desired accuracy. (Default: tol = 1.e-6).
// trace - If nonzero, each step is printed. (Default: trace = 0).
//
// OUTPUT:
// out - Returned integration time points (1st column).
// Remaining columns are solutions of equations specified in Ufun
//
// The result can be displayed by: plot(tout, yout).
// Daljeet Singh
// Dept. Of Electrical Engg., The University of Alabama.
// 11-24-1988.
// Modified by Ian Searle
// 4-18-93
static (collapse);
ode78 = function (Ufun, t0, tfinal, y0, Uoutput, tol, trace)
{
local (alpha, beta, chi, delta, f, gamma1, h, hmax, ...
hmin, j, psi, pow, t, tau, tout, y, yout, lout, I, output);
// The Fehlberg coefficients:
alpha = [ 2./27., 1/9, 1/6, 5/12, .5, 5/6, 1/6, 2/3, 1/3, 1, 0, 1 ]';
beta = [ 2/27, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
1/36, 1/12, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
1/24, 0, 1/8, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
5/12, 0,-25/16, 25/16, 0, 0, 0, 0, 0, 0, 0, 0, 0 ;
.05, 0, 0, .25, .2, 0, 0, 0, 0, 0, 0, 0, 0 ;
-25/108, 0, 0, 125/108, -65/27, 125/54, 0, 0, 0, 0, 0, 0, 0 ;
31/300, 0, 0, 0, 61/225, -2/9, 13/900, 0, 0, 0, 0, 0, 0 ;
2, 0, 0, -53/6, 704/45, -107/9, 67/90, 3, 0, 0, 0, 0, 0 ;
-91/108, 0, 0, 23/108,-976/135, 311/54, -19/60, 17/6, -1/12, 0, 0, 0, 0 ;
2383/4100, 0, 0, -341/164, 4496/1025, -301/82, 2133/4100, 45/82, 45/164, 18/41, 0,...
0, 0 ;
3/205, 0, 0, 0, 0, -6/41, -3/205, -3/41, 3/41, 6/41, 0, 0, 0 ;
-1777/4100, 0, 0, -341/164, 4496/1025, -289/82, 2193/4100, ...
51/82, 33/164, 12/41, 0, 1, 0 ]';
chi = [0, 0, 0, 0, 0, 34/105, 9/35, 9/35, 9/280, 9/280, 0, 41/840, 41/840]';
psi = [1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, -1, -1 ]';
pow = 1/8;
if (!exist (tol)) { tol = 1.e-6; }
if (!exist (trace)) { trace = 0; }
// Initialization
t = t0;
// the following step parameters are used in ODE45
// hmax = (tfinal - t)/5;
// hmin = (tfinal - t)/20000;
// h = (tfinal - t)/100;
// The following parameters were taken because the integrator has
// higher order than ODE45. This choice is somewhat subjective.
hmax = (tfinal - t)/2.5;
hmin = (tfinal - t)/10000;
h = (tfinal - t)/50;
y = y0[:];
f = y*zeros (1,13);
tout = t;
yout = y.';
if (!exist (Uoutput))
{
output = 1;
lout = << [tout, yout] >>;
else
output = 0;
lout = << Uoutput (t, y).' >>;
}
I = 2;
tau = tol * max( [norm(y, "I"), 1]);
if (trace) { [t, y[:].'] ? }
// The main loop
while ((t < tfinal) && (h >= hmin))
{
if (t + h > tfinal) { h = tfinal - t; }
// Compute the slopes
f[;1] = Ufun (t, y)[:];
for (j in 1:12)
{
f[;j+1] = Ufun (t+alpha[j]*h, y+h*f*beta[;j])[:];
}
// Truncation error term
gamma1 = h*41/840*f*psi;
// Estimate the error and the acceptable error
delta = norm (gamma1,"I");
tau = tol*max ([norm (y,"I"),1.0]);
// Update the solution only if the error is acceptable
if (delta <= tau)
{
t = t + h;
y = y + h*f*chi;
if (output)
{
lout.[I] = [t, y.'];
else
lout.[I] = Uoutput (t, y).';
}
I++;
}
if (trace) { [t, y[:].']; }
// Update the step size
if (delta != 0.0)
{
h = min ([hmax, 0.8*h*(tau/delta).^pow]);
}
}
if (t < tfinal)
{
fprintf ("stderr", "SINGULARITY LIKELY, at t = %d\n", t);
}
return collapse (lout);
};
//
// Collapse a LIST of matrices into a single matrix.
//
collapse = function ( LIST )
{
local (i, m, row, col);
row = size (LIST.[members (LIST)[1]])[1];
col = size (LIST.[members (LIST)[1]])[2];
m = zeros (length (LIST)*row, col);
for (i in 1:length (LIST))
{
m[i;] = LIST.[i];
}
return m;
};